from rmpoly import Poly

def det_rec(m):
  """computation of the determinant in a recursive way
  reducing it to triangular form
  taken from B.M.Bingle 'Symbolic determinants:calculating the degree'
  >>> from rmpoly import *
  >>> from fractions import Fraction as mpq
  >>> print det_rec([[1,2,3],[mpq(1,2),0,4],[0,2,4]])
  -9
  >>> rp,y,x = rgens('y,x',6,mpq)
  >>> det_rec([[x+y, x],[y,x-y]])
   +x^2 -x*y -y^2
  """
  n = len(m)
  if n == 1:
    return m[0][0]
  if n == 2:
    return m[0][0]*m[1][1] - m[0][1]*m[1][0]
  # construct submatrix (with 1 less row and column than m)
  # if first row starts with 0, exchange it with the first
  # row not starting with 0
  i = 0
  while 1:
    if m[i][0]:
      break
    i += 1
  if i == n:
    return 0
  if i:
    m[0],m[i] = m[i],m[0]
  sign = (-1)**i
  m1 = []
  for i in range(1,n):
    a = [m[i][j]*m[0][0] - m[i][0]*m[0][j] for j in range(1,n)]
    m1.append(a)
  det = det_rec(m1)
  den = m[0][0]**(n-2)
  if sign == 1:
    return det/den
  else:
    return -det/den

det = det_rec

def neville_interpolation(x,points):
  """univariate polynomial interpolation with Neville algorithm
  x variable; x.rp.bits_exp must have enough bits for the polynomial
  points = ((x0,y0),(x1,y1),...) such that p(x=x_i) = y_i
  see http://en.wikipedia.org/wiki/Neville%27s_algorithm
  >>> from rmpoly import *
  >>> from fractions import Fraction as mpq
  >>> rp,x = rgens('x',6,mpq)
  >>> neville_interpolation(x,((0,2),(1,5),(2,10)))
   +x^2 +2*x +2
  """
  n = len(points)
  field = x.rp.field
  xv = [field(points[i][0]) for i in range(n)]
  yv = [field(points[i][1]) for i in range(n)]
  m = [[0]*n for i in range(n)]
  for i in range(n):
    m[i][i] = yv[i]
  for k in range(1,n):
    for i in range(n-k):
      j = i+k
      m[i][j] = ((x-xv[j])*m[i][j-1]+(xv[i]-x)*m[i+1][j])/(xv[i]-xv[j])
  return m[0][n-1]




if __name__ == "__main__":
    import doctest
    import sys
    if sys.version_info < (2, 6):
      print 'doctests require Fraction, available from Python2.6'
      sys.exit()
    doctest.testmod()

